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We generalize the momentum average (MA) approximation to study the properties of models with 
Pq momentum-dependent electron-phonon coupling. As in the case of the application of the original 
MA to the Holstein model, the results are analytical, numerically trivial to evaluate, exact for both 
zero bandwidth and for zero electron-phonon coupling, and are accurate everywhere in parameter 
space. Comparison with available numerical data confirms this accuracy. We then show that 
further improvements can be obtained based on variational considerations, using the one-dimensional 
breathing-mode Hamiltonian as a specific example. For example, by using this variational MA, we 
obtain ground state energies within at most 0.3% error of the numerical data. 
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I. INTRODUCTION 

There is considerable need to understand the coupling 
of a particle to its environment, in particular the interac- 
tion between a charge carrier and phonons. The physics 
of polarons - electrons dressed by phonon clouds - is be- 
lieved to be relevant in explaining a plethora of phys- 
ical systems and properties, for exa mple in polymers, 
nanotubes. Ceo and other fullerenes,'32ll! manganites,* 
Bechgaard salts, 1^ and even, possibly, the kink in the 
quasiparticle dispersion in cuprates, a necessary ingre- 
dient for the ful l unde rstanding of high-temperature 
superconductivity.'^'^'Slini 

It is of obvious benefit to have highly accurate ana- 
lytical approximations that can be applied to compli- 
cated many-body problems, such as the polaron prob- 
lem. Unfortunately, most existing approximations fail to 
reproduce the correct physics in "intermediate" regimes 
wher e the relevant physics often occurs. In recent 
workd^3iIIilElIiIIIll we have proposed the so-called mo- 
mentum average (MA) approximation for calculating the 
Green's function of a dressed particle, focusing on the 
case of an electron dressed by phonons. This approx- 
imation is analytical, easy-to-use, and highly accurate 
throughout all of parameter space, both at low and at 
high-energies. Its underlying idea is to sum all of the 
diagrams in the diagrammatical expansion of the self- 
energy, but with the diagrams approximated in such a 
way that the full summation can be performed. At the 
simplest level, this is achieved by replacing all of the free 
propagators appearing in the self-energy by their momen- 
tum average. For the Holstein model, the resulting MA 
self-energy exhibits not only the exact asymptotic be- 
haviour for zero coupling and zero bandwidth, but it is 
also in excellent agreement with numerical results in the 
intermediate r egimes w here other approximate methods 
completely fail.^HEIl 

In Refs. [TT] and [T^ we presented a detailed account 
of the derivation and application of the simplest-level 
MA approximation (now known as MA*^°-') to the Hol- 
stein model, including extensive comparisons to the avail- 



able numerical data. We further justified the accuracy 
of the approximation by looking at the spectral weight 
sum rules. The MA^'^^ approximation satisfies the first 
six spectral weight sum rules exactly, and is highly accu- 
rate for all higher order sum rules. In Ref. [T5]we also 
pointed out the three key limitations of the MA'^'^^ ap- 
proximation: (i) it fails to correctly predict the so-called 
electron+phonon continuum that must occur at a phonon 
energy above the ground state energy of the polaron, (ii) 
its accuracy worsens for small phonon energies, and (iii) 
the MA^''^ self-energy is independent of momentum. 

These shortcomings motivated the systematic improve- 
ment of the MA approximation, presented in Ref. [TH 
This lead to a hierarchy of increasingly more accurate ap- 
proximations, which we call MA(°), MA^^^, MA^^^ etc. 
The increased accuracy is a result of fewer approxima- 
tions for the self-energy diagrams, however done in such 
a way that we can still sum all the resulting diagrams 
analytically. It is also possible to understand the MA ap- 
proximations in a variational context, as detailed below. 
From this point of view, the systematic improvements 
are obtained by including additional states required to 
reproduce correctly the electron+phonon continuum. In 
the process we also obtain a momentum-dependent self- 
energy and higher accuracy throughout all of parameter 
space, particularly for small phonon energies. As a re- 
sult, the MA approximations allow us to understand very 
accurately the Holstein polaron physics, throughout the 
parameter space, for all energies and momenta. 

In this paper we generalize this powerful set of approxi- 
mations to a much broader class of models with electron- 
phonon coupling that depends on the phonon momen- 
tum. We first derive a simple generalization, leading to 
what we will continue to call the MA(°), MA^^), MA^^), 
etc. hierarchy. These are very easy to apply to any 
Hamiltonian of this class, however while still asymptoti- 
cally exact for both weak and strong coupling, at inter- 
mediary couplings the relative errors are of a few percent. 
In other words, these can be used to get a quick estimate 
of typical energies and spectra. We then show how these 
can be significantly improved, using variational consider- 
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ations. This generates a second hierarchy of approxima- 
tions which we call MA^''-*", MA^^^^), MM'''^\ etc., with 
relative errors well below one percent. As a test case 
to gauge these accuracies, we use the one-dimensional 
(ID) breathing-mode Hamiltonian, where high-accuracy 
numerical results have recently become available.!^ 

The general electron-phonon coupling model that we 
consider has the following form in momentum space: 
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The first term is the kinetic energy of the electron, where 
cj^ and Ck are electron creation and annihilation opera- 
tors, and Ek is the electron dispersion. For the single 
electron (polaron) problem of interest to us, the spin 
of the electron is irrelevant and we suppress its index. 
The second term describes a branch of optical phonons 
of energy f2, 6^^ and bq being the phonon creation and 
annihilation operators. The last term describes the cou- 
pling of the electron to the phonons, where gq is the 
momentum-dependent coupling. Sums are over all mo- 
menta inside the first Brillouin zone, and we set h = I 
and a — 1 throughout this paper. This general Hamilto- 
nian covers complicated electron-phonon couplin gs s uch 
as those found i n the Rashba-Pekkar,'l^li2l Frohlichp^l and 
breathing-mod^^^^'^D models, and it reduces to the Hol- 
stein model when is simply a constant .122! 

The quantity of intere st to us is the (retarded) single 
polaron Green's function! ^^ * ^^ ! 
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where |0) is the vacuum, Ck|0) — bq\0) — 0, and 77 > 
is infinitesimally small. The importance of this Green's 
function is obvious in the Lehmann representation;^'^ 
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full 2D form, it describes lattice vibrations in a CUO2- 
like plane, where the motion of the O ions living on the 
bonds connecting the Cu sites is the most important vi- 
brational degree of freedom,'^ making it possibly rele- 
vant for the study of high-Tc superconductors. The sec- 
ond reason is that exact diagonalization (ED) result^^^ 
have recently become available for its ID analog, relevant 
for CuO chains. We focus here on this ID breathing- 
mode model because these results serve as an excellent 
gauge of the accuracy of the generalized MA approxi- 
mations. Their availability is very fortunate because al- 
though there are many numerical results for the Holstein 
model, it is only due to recent advancements in computa- 
tional power that more complicated electron-phonon cou- 
pling models, such as the breathing-mode Hamiltonian, 
can be investigated numerically. We also mention that 
this model has been studied using the self-consistent Born 
approximation (SCBA),'2iIbut SCBA is known to b e very 
poor for intermediate and large coupling strengths.l^^EZl 

In the ID breathing-mode model, one considers a chain 
with two interlaced sublattices, where the Cu sites which 
host the electron are indexed by integer labels, and the 
O sites which host the phonons are indexed by half- 
integer labels. The interaction term of the breathing- 
mode Hamiltonian can be written in real space as 
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1^2) where Xi±i/2 = ^i±i/2+^i±i/2 describe the displacements 



of the O atoms neighboring the Copper atom at site i 
that hosts the electron, and g is a constant describing 
the strength of the electron-phonon coupling (g absorbs 
the proportionality factors between the true displacement 
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2 and the phonon operator b 
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). Trans- 



forming into momentum space, the ID breathing- mode 
Hamiltonian takes the form of Eq. ([ij with 



where {|a)} and {E^} are the complete set of one-particle 
eigenstates and eigenenergies, Tl\a) = Ea\a). In this rep- 
resentation it is clear that the poles of the Green's func- 
tion give the one-particle spectrum, and the associated 
residues, sometimes called quasiparticle (qp) weights, 
give partial information on the nature of the eigenstates. 
Knowledge of this Green's function also allows us to cal- 
culate other relevant quantities such as the effective mass 
of the polaron, the average number of phonons in the 
polaronic cloud,'^^^!! or more detailed phonon numbers 
statistics.'^S] Furthermore, its imaginary part is directly 
measured experimentally by angle-resolved photoemis- 
sion spectroscopy (ARPES).— 

In this paper we will use the ID breathing-mode Hamil- 
tonian as an example, however, we stress that the meth- 
ods we present are applicable to any Hamiltonian like 
Eq. in any dimension. The reason for choosing the 
ID breathing-mode Hamiltonian is two-fold. First, in its 
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In this model the electron motion is described by a tight- 
binding model with the usual e^, = —2t cos{ka) , although 
our results can be applied for any dispersion. 

This article is organized as follows. In Sec. II, we de- 
rive the exact set of equations giving the Green's function 
of a polaron, review the MA approximations and derive 
the simple straightforward generalizations MA*^°\ MA^^^ 
etc. for models similar to Eq. ([ij. In Sec. HI we show 
how to obtain the more accurate MA^""'*^) and MA^"''^) 
approximations based on variational ideas. Here we use 
the ID breathing-mode Hamiltonian as an explicit exam- 
ple. In Sec. IV we present our results and compare them 
to the available numerical data where possible. Finally, 
Sec. V contains our summary and conclusions. 
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II. CALCULATING THE GREEN'S FUNCTION 

A. Exact Solution 

Although exact solutions for these types of problems 
are generally not obtainable in closed form, one can 
write down their formal solutions in terms of an infinite 
set of coupled equations involving related (higher-order) 
Green's functions. In a previous worlf^^ we described in 
detail how to generate these equations for the Holstein 
model. The generalization to momentum-dependent 
coupling models is straightforward, by repeated use of 
Dyson's identity G{uj) = Goi^) + G{lj)VGq{uj) , where 



G^[uj~'H + ir]]-\ Go = [Lu-no + ir)]-\'H = 'Ho + V, 
and V is the electron-phonon interaction. Following this 
procedure and defining the generalized Green's functions 

F„(k, qi, . . . , q„, CO) = {OKG{iu)cl_^X^ ' ■ ' l^) (6) 
where Fo(k,[x)) = G(k,cj), we find that 



G(k,c^) = Go(k,c^) 



and for ?i > 1, 
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V A* 



.1=1 

+ ^ 9cii+iFn+i{k,qi,...,qn+i,i^) 



qn+i 



(8) 



The total momentum carried by the phonons is denoted 
by qT = EILi q» a'^d Go(k, w) = (w - et + is the 

free electron Green's function. Observing from Eqs. ([s]) 
that all of these generalized Green's functions Fi , F2 , . . . 
must be proportional to G(k,a;),'^ we can recast our 
equations into a more convenient form by defining 



/«(qi 



^ N _ ^"^^ffqi ■■■gq„-F„(qi,...,q„) 
G(k,c^) ' 



where we have also introduced the shorthand notation 
/„(k,qi, . . . ,q„,Lj) = /„(qi, . . . ,q„) (i.e., the k and uj 
dependence of these functions is implicitly assumed from 
now on). In this notation, Eq. ([t]) becomes 



G(k,c^) = Go(k,c^) 



i + ^E/i(qi)G(k,c.) 



with a solution written in the standard form 
G(k,w)= ^ 



Lu — e^ — S(k, Lu) - 
where the self-energy is 

S(k,a;)=l^/i(qi). 
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This self-energy can also be written in terms of an infinite 
set of diagrams, as shown in Fig. [T] 

The quantity of interest /i(qi) is obtained from the 
set of equations resulting from Eq. (pi , namely /o = 1 , 



by definition, and for n > I: 

/«(qi,---,q«) = Go{k- qT,uj - nfl) 

71 

^ l5q.P/"-i(- ■ • >q«-i>q»+i: • ■ •) 
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Of course this system can be solved trivially in the 
limit of g = 0, in which case G(k, w) = Go(k, w) directly 
from Eq. ([7|. Also, in the limit of t = the free propa- 
gators become independent of momentum and the equiv- 
alent of the Lang-Firsov result is reproduced.'^ However, 
for the general case of finite g and t, a closed form solu- 
tion cannot be obtained, even for the Holstein model.l^ 
Approximations are therefore needed. We begin by de- 
scribing the simplest MA^^^^ version of a generalized MA 
approximation. 
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FIG. 1: Diagrammatical expansion of the self-energy S(k, c 
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B. The MA^°' Approximation 

For the Holstein model, the MA^^^ approximation 
amounts to replacing all of the free electron propaga- 
tors in the diagrammatical expansion of the self-energy 
by their momentum average over the Brillouin zone: 



(14) 



which is equivalent to replacing all G'o(k — qr,^ — nfl) 
by goiio — nfi) in Eq. (13). The procedure is essentially 



the same for momentum-dependent electron-phonon cou- 
plings. We note that the first term from the n = 1 case 
of Eq. ( 13 ) does not actually require any approximation 



(this is trivially true in the case in the Holstein model) 
because /o = 1 by definition, and its coefficient can be 
written explicitly as 



ffo(k,a;) = ^^|5q|'Go(k-q,w). 



(15) 



By making the substitution Go(q, — nCl) — > gQ(LU — nfl) 
everywhere else, we can rewrite Eqs. (12 1 and (IT3|) in 



terms of the momentum averaged functions JF„(k,ti;) — 
h Sqi,...,q„ ■ • ■ ' q™)- This leads to simple recur- 
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rence relations linking each JF„ to Tn-i and Tn+i, and 
these can be solved in terms of continued fractionsJ ^^ * ^*^ ! 
We find 



(k,^) 
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where we have defined the infinite continued fractions 
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and the momentum averaged electron-phonon coupling: 

1 



N ^ 
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It is trivial to check that this reduces to the Holstein 
result of Ref. [T^when g^ is independent of momentum. 

Before we discuss how to systematically improve this 
result, we briefly review one explanation as to why this is 
a reasonable first step in obtaining an accurate approxi- 
mation: in real space this approximation is equivalent to 
replacing, in all self-energy diagrams, all free propagators 
Go (i , j , w — n fl) — > SijgQ{uj — nfl), where i and j index the 
electron sites.'^^HH jf ^j^g fj-gg propagators Go{i,j,uj — n^l) 
are evaluated for energies well below the free electron 
continuum (and for a system with interactions the po- 
laron ground state is below the free electron continuum). 



it is well known that the free propagator decreases ex- 
ponentially with increasing distance \i — j\. Thus, the 
most important terms are those with i = j, i.e. pre- 
cisely those included within MA*^"), explaining why the 
approximation should be reasonably accurate at low en- 
ergies (spectral weight sum rules then insure that it is 
similarly accurate for all energies). 

This also suggests a means to systematically improve 
the MA^'^^ approximation. Since the propagators with 
the higher energy (the energy closer to the free-electron 
continuum) decrease exponentially with distance most 
slowly, the first improvement is to keep all propagators 
with energy to — exactly. Higher order systematic im- 
provements denoted by MA^") would amount to keeping 
all propagators with the argument to — mfl exactly, for 
all m < n. In Ref. [THwe derived explicitly the Holstein 
MA self-energy for both n = 1 and n = 2. As shown 
and explained there, the n = 1 order already insures the 
key improvement of properly predicting the polaron-|-one 
phonon continuum, which is usually absent at the n — 
level. The reason we went to n — 2 for Holstein is that, 
given the simplicity of that model, only at the MA^^^^ 
level did we find a momentum-dependent self-energy. 

In contrast, for a model with a momentum-dependent 
coupling even the MA(°) level gives explicit k dependence 
in the self-energy [see Eq. (16l], therefore in the next 
subsection we only consider the MA'^^ generalization. 



C. The MA^^' Approximation 



The main drawback of the MA^°^ level of approxima- 
tion is that it predicts an incorrect location for the elec- 
tron+phonon continuum that must start at Eqs + f2. 
This problem is always "cured" at the MA^-'^^ level be- 
cause in real space MA^^^ includes states with a phonon 
cloud near the electron and a single phonon arbitrarily 
far away, in other words precisely the type of states that 
give rise to the polaron-|-phonon continuum.^ Mathe- 
matically, MA^^^ amounts to keeping all propagators with 
the argument lu — fl in Eq. ( 13 1 exactly. Such terms only 



appear in the n = 1 equation: 



/i''^(qi) = Go(k-qi,u;-f)) \g^,\ 



(qi,q2) 

(19) 

where we have distinguished the approximated /„ terms 
with the superscript (1) to indicate the MA^^-' level of 
approximation. For the remaining equations (n > 2) we 
proceed as before, replacing Go(k — qx,'-^ — nfl) with 
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(7o(w — nf2) everywhere in Eq. ( 13 1: 
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(20) 



We wish to solve for EJ^^,^,{k,iu) = (l/N) J2^^ f['\m)- 
The procedure is analogous to that of Ref. [TH We 
obtain two sets of coupled recurrence relations, one 
for fully momentum averaged quantities ^„(k,[j) — 

.,q„) that already appeared at 
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MA(o) level, and one for partially averaged quantities 

:r„(qi,k,w) = ]v^^Eq2,...,q,. /«^^(qi>--->qn)- Their 
solution follows the procedure detailed in Ref. [TH and 
we simply state the result: 



^A/A(i)lK,CJj — ^— — ^ — 

1 - 5o(k,w) [A2(uj) - Ai{uj - 

where uj = u—fl—g^Ai^uj—fl). This expression is slightly 
more complicated than the expression for Ejyj^co) (k, uj) 
because it involves two continued fractions, but it is again 
trivial to compute numerically. 



D. MA(°) and MA^^^ Results 

To illustrate the accuracy of MA(o) and MA^^) we will 
use the ID breathing-mode Hamiltonian described by 
Eqs. ^ and (|5| as an example, comparing our results 
to those found numerically using ED in Ref. |T71 

Before reporting these results we briefly discuss how 
the MA results are evaluated in practice. First, the nu- 
merical evaluation of the infinite continued fraction found 



in Eq. ( 16 1 requires that we truncate the fraction at some 
level n. For an error of order e we require that at the n^^ 
level we have e > ng^go^uj — nfl)gQ{uj — (n + 1)51). Since 
go{uj) — > l/w when |w| — > oo, for large enough n we can 
approximate 50(1^ ~ nfl) ~ .go('^ — {n + 1)51) ~ —l/(nfl). 
It then follows that for an error of order e, we must have 



n > 



(22) 



In practice we always check our results by doubling the 
value of the truncation value n until the change in the 
self-energy is negligible. The result is that all MA error 
bars throughout this paper are smaller than the thickness 
of the lines/symbols used in the plots. 

With our approximations for the Green's function in 
hand, the ground-state properties are found by tracking 
the energy and weight of the lowest pole of the spectral 




-4.5 

5 

g4 

~8 3 

0: 





_ ■ — ■ ma'" 

- •--•ma'"' 

- ^< 


, 1 1,4, 1 1 1 , 

/ (b) - 

* J ^^-^ 




JtT 1 , , , 1 , , , " 



0.5 1 1.5 

A,=g^/(at), X=X/2 



0.5 1 1.5 

A.=g'(nt) 



FIG. 2: (Color online) (a) Ground state energy, (b) Percent 
difference from ED ground state energy results, (c) Ground 
state qp weight, and (d) Percent difference from ED qp 
weight results, as a function of the effective coupling A, for 
t = l,Q/t — 0.5. The perturbation theory results^- for both 
the Holstein and breathing-mode models are shown, and the 
Holstein result is plotted as a function of Ah = A/2. The ED 
results are from Ref. 1171 



weight, defined as 
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(23) 

where a small but finite value of rj broadens the 6 peaks 
into Lorentzians and enables us to detect them numeri- 
cally (typically we take rj ~ 10~^). A detailed descrip- 
tion on how we use the Green's function to extract the 
energy spectrum and quasiparticle (qp) weights, as well 
as other quantities of interest such as effective masses, 
average number of phonons in the polaron cloud, etc., is 
presented in Ref. 1121 The explicit expressions for go{uj) 
and gi){k,Lo) needed to compute the self-energies for the 
ID breathing-mode Hamiltonian are given in Appendix 
\K\ As is customary, we define the dimensionless effective 
coupling A as the ratio between the lattice deformation 
energy ^7 —2g'^/^ and the free electron ground state en- 
ergy -2t: 



nt' 



(24) 



In Fig. [2] we plot the ground-state energy and qp 
weight as a function of the electron-phonon coupling 
strength, for a phonon energy fl/t — 0.5. The MA'^'^^ (red 
dashed line) and MA^^^ (solid blue line) results both show 
good agreement with the ED results (black circles). The 
ground state energies calculated with MA*^") are within 
5% error of the ED results, and the MA^^^ results are bet- 
ter, coming within 3% error of the exact energies. They 
are exact for both t = and A = 0, as expected, and 
the crossover from the large to small polaron is captured 
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to a high degree of accuracy. This accuracy is very en- 
couraging, especially since these approximations are so 
trivial to evaluate. It is also worth pointing out that 
our work on the Holstcin Hamiltonian shows that this 
accuracy improves in highcr-dimcnsional models, and we 
believe this to be true here as well. Unfortunately, lack of 
detailed numerical results in higher dimensions prevents 
us from confirming this to be the case for models with 
momentum-dependent coupling. 

Two more observations are apparent regarding these 
results: (i) The energies predicted by MA are lower in 
energy than the ED results, and (ii) the MA results ap- 
proach the Lang-Firsov asymptotic limit very slowly. 

The fact that the energies predicted are below the ex- 
act result indicates that the MA approximation is non- 
variational in the case of the breathing-mode model. This 
is somewhat surprising, since for the Ho lstein model it 
has been shown that MA'") is variationaP^EH and that 
MA^^' is quasi-variational (in the latter case the Hamil- 
tonian is modified slightly and the approximation is no 
longer truly variational, for details see Ref. [H]) . In any 
case, the energies found for the Holstein model using the 
semi-variational MA*^^^ and MA*^^) approximations were 
always higher than the exact numerical results. 

The slow asymptotic convergence at large A is due to a 
difi^erent prefactor of the 0{P) perturbational correction. 
Instead of the correct breathing-mode result!^ 



2tcos{ka)e 
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(25) 



we have produced the Holstein result from Ref. [T2l with 
replaced with g^ = 2g^ [see the two dashed lines la- 
beled PT in Fig. ^a.)]: 
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(26) 
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To understand the origin for both of these facts, we 
observe that at this level of approximations, almost 
all dependence on the el-ph coupling (7q is through its 
momentum-averaged value cf^ — l/A^X]qlffgP- For the 
ID breathing-mode model, this average happens to be 

the same whether cx isin^, as is the case here, or ([l6| and (|2T]) are very easy-to-use, and rather accurate 
whether gq oc cos ^ which would correspond to a cou- 



FIG. 3; MA*"' and MA'^' predictions of the ground-state en- 
ergy of the HIC modeP^ as compared to its exact solution 
(black dots). For this model with infinite-range electron- 
phonon coupling, the approximations capture the correct 
asymptotic behavior. 



towards the charge. The effective hopping is related to 
the overlap of these clouds when the charge hops to a 
neighboring site. If the electron hops from i to i + 1, the 
phonon cloud at i -|- ^ changes from exp ^&|^i]|0) to 

exp 1 ] |0) and the overlap is small. In the MA ap- 
proximations, it is precisely this information about left 
vs. right that is lost (equivalently, one can think of los- 
ing the information about relative phases between vari- 
ous contributions to the phonon cloud) and the overlap 
is unity - both cases have a cloud on the i + \ This 
explains the higher polaron mobility, hence the lower en- 
ergy and slower convergence towards the exact asymp- 
totic value for the MA approximations. 

For the Holstein model, this problem does not appear 
because the el-ph coupling is local. We argue that this 
problem should also become less serious if the el-ph cou- 
pling is longer range, because in that case the relative 
sign of various displacements is not completely lost in 
the — average. Unfortunately, lack of 

numerical data for lattice models with longer-range el- 
ph coupling makes it difficult to support this statement. 
The only comparison we can offer is with the exact solu- 
tion of a somewhat pathological, infinite-range coupling 
model called the HIC model of Ref. [27l As shown in Fig. 
[s] for this model both MA*^"^ and MA^-'^^ approximations 
give the correct asymptotic behavior - and this necessi- 
tates a correct description of the infinite-range phonon 
clouds that appear in this model. 

To summarize, the MA approximations given by Eqs. 



pling of the electron to the sum x^^^i + x^j^i of 0-site 
displacements. In other words, for this model, these 
simplest versions of the MA approximation register that 
a phonon cloud is formed at a certain O site, but not 
whether this leads to a leftwards or rightwards displace- 
ment of that site. In the strong coupling limit, one ex- 
pects clouds to be formed only on the two O sites neigh- 
boring the Cu site that hosts the charge, and to point 



for models with momentum-dependent electron-phonon 
coupling. Given its low dimension and short-range (but 
not local) el-ph interaction, one would expect the ID 
breathing- mode model to be amongst the worst examples 
for the accuracy of this approximation. However, even 
here we obtain very decent agreement. As argued, we 
expect this to improve for higher dimensions and longer- 
range interactions. The only regime where accuracy is 
certain to worsen is when f2/i becomes very small (this is 
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FIG. 4: The generalized Green's functions that appear in 
ma'"' are sketched in the top picture. They have phonon 
clouds on two neighboring sites, with the electron on the 
central site. For the variational MA approximation, denoted 
we find again such Green's functions, but also those 
with the electron to the immediate left or right of the phonon 
clouds. We name them f^m and fnmj respectively. 



a general problem of MA- like approximations, discussed 
at length in Ref. [T^ . Therefore, we believe that these 
approximations are useful for quick guides to relevant 
energies and other quantities of interest for problems of 
this type, obtained with minimal analytical and compu- 
tational effort, yet still reasonably accurate. 



III. THE VARIATIONAL MOMENTUM 
AVERAGE APPROXIMATION 

In this section, we attempt to remedy the problems 
pointed out in the previous section and by so doing, 
to obtain an improved MA approximation for the ID 
breathing-mode Hamiltonian. The solution, which is pre- 
sented in this section, can then be used as a template for 
any other model. The main idea is to try to formulate 
an MA approximation which is variational in n ature, as 
the original MA is for the Holstein Hamiltonian.fSHH 



A. The MA^"'"' Approximation 



a guidance for what states are most relevant to in- 
clude. Remember that this solution involved only 
the fully momentum-average quantities J-"„(k, lj) = 
Sq, ,...,q,. /n(k, qi, . . . , q„, w). Using the definition of 
Eqs. (|9j and ([6| and performing the sums over gi, . . . q„ 
for our specific model with gq given by Eq. (|5|, we 
find immediately that J2i e*'^^' (0|ckG'(u;)c|(6|_ i — 

foj_i_i)^|0). In other words, only states where there are 
phonons only on two neighboring sites contribute to it. 

We use this as the criterion for our variational MA ap- 
proximation. More specifically, when we use the Dyson 
equation to generate equations of motion, we only keep 
the terms which have phonons only on two neighboring 
sites, and discard any other contribution. Of course, 
these states have to have a total momentum fc, or else 
the matrix element is zero. We find that 3 sets of gener- 
alized Green's functions, sketched in Fig. |4] appear when 
we use this criterion, namely: 



NG{k,Lj) ^ \ I V ; ^ ,_i ,+ 1 I 



(27) 



and 



7VG(fc,tj) 

(28) 

Note: as before, we again we do not write explicitly the 
fc and Lo dependence. 

It is straightforward to show that in this notation we 
can write Eq. ( 10 1 as 



G(fe,t^) = Go(fc,^)[l -H5(/i,o - /i,i)G(fc,a;)], (29) 



and therefore 



(30) 



Given the good agreement found in the previous sec- 
tion, it is reasonable to use the MA^^^ solution as 



The equations of motion consistent with the variational 
restriction are straightforward to find: 



/n,0 — eo \nfn-\Si + /n+1,0 " /n+1, 



626 



J n-t-l,n 

— 6i \ne~ 



Jn+l,n+l 



(31) 
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For any 1 < m < n — I, 

fn,m — Co [in — m)fn-l,m — mfn-l^rn-l + fn+l,m — /n+l,m+l] 

— ei [(n — 'Ti)/^_j^ „j — m/„_j^_„j__]^ + fn+l,m ^ /ri+l,m+l] ' (32) 

and 

fn,n = —Co [f^/n_l,n-l + fn+l,n+l — fn+l,n] + ^-2^-^''°' .fn+1,1 

+ ei [ne'^^fn-ifl - e*""/„+ia + e^'Vn+i.o - /++i,„] • (33) 

Similarly, for any 1 < m < n ~ 1, 

fn,m — [('^ ~ rnjfn-l^m ~ ^fn-l,m-l + fn+l.m ~ fn+l.m+l] 



eo 

62 



[in - m)f+_^ ,^ + /++i^„J + I eo } [™/«-i,m-i + /«+i,m+i] (34) 



and 

fn,0 = 6i [nfn-lfi + /n+1,0 - /n+l,l] + '^ofn+1,1 ~ ^3^^^'"' fn+l,n 

- 62 [ne /„-i,„-i + e /„+i,n+i - e /n+i,„J (35) 



'61 [n/„„i,„_i + /„+!,„+! - /n+l,n] - 6o/„+i_„ + 636* °/„+i_i 



62 [ne'"''!, 



71-1.0 " 



yn+1,0 



+14 



(36) 



where we use the shorthand notation ej = ggj{ijJ — n^l), 
and gj{uj) are real-space Green's functions, defined as 



9j 



H ^ ^^6±»«(^'^)Go(g,u;) = Go(j,0,^). (37) 

q 



The "±" sign in the exponent of Eq. (37 1 is irrele- 
vant because GQ{q,Lo) is even with respect to q. The 
explicit expressions for these and other momentum aver- 
aged functions of the free electron propagator are given 
in Appendix [A| for a tight- binding dispersion. Note that 



/n,n = 6'''V«,o and 



/n,o, which is why we do 



not need to keep them as independent variables. 



Equations (31 )-(36l can now be cast in the form v„ = 
+ B„v„_(-i, where 

Vn = (/n.Oj • ■ • I /n,n: fn.O^ ■ ■ ■ ' /n.n-li • ■ • ' /n^n) 

(38) 

collects all generalized Green's functions with a total of n 
phonons, and the matrices A„,B„ are straightforward to 
obtain from the equations above. The solution of this set 
of recursive equations can then be written as an infinite 
continued fraction involving products of matrices?^ 



1 



1 -B 



-A„v„_i. 



(39) 



Although this is a continued fraction of matrices of in- 
creasing size, it is still numerically trivial to evaluate. 



The dimensions of A„ and B„ are (3n-|- 1) x (3n — 2) and 
(3n -I- 1) X (3n -t- 4), respectively, and we determine the 
truncation level of the continued fraction using the same 
criteria as discussed in the previous section for the MAf^") 
approximation. In addition, the A„ and B„ matrices are 
very sparse, which makes multiplication by them very ef- 
ficient. Finally, since wq = 1 by definition, it is straight- 
forward to calculate Ej;,j^(^,o) (A:, w) from Eq. (30 1 once 
fiO'fti)'^ = Qi is evaluated from Eq. 



Vl 



(/i,o, /i 



(391 



B. MA'"'^' approximation 



We can also systematically improve the variational MA 
approximation to reproduce the electron-Kphonon con- 
tinuum. It is obvious that this will not be predicted by 
since no phonons are allowed to appear far from 
the main polaronic cloud. 

We follow the same approach as before: at the MA^^^'^^ 
level we keep all equations involving free electron prop- 
agators with u! — U exactly. In order to work in the 
enlarged variational space described by Eqs. (27) and 



(28 1 we need to define the following generalized Green's 



functions (these are analogous to the partial momentum 



9 



averages ^(qi,k,tj)): 



/n,m('Zl) 

and 

fn,mi.<ll) 



9qi 



,j(fc-9i)i?.i 



7VG(fc,w) 

X (0|c,G(c.)4&t^6t™ fe|(T''"'"|0) (40) 



^i(k-qi)Ri 



VNG{k,uj) 
X (0|c,G(c.)c,Li6Ut™ ^t(-i)-™|o). (41) 



These Green's functions explicitly contain states with one 
phonon delocalized away from the main polaronic cloud, 
i.e. precisely the type of states required to reproduce the 
polaron+phonon continuum. In this notation Eq. ( 19 ) 
can be written as 



/i(<7i) - Go{k~q,,LU-n) [\g,f + 5(/2.o(gi) - hAn))] . 

which is an exact equation involving no approximations, 
as in the MA'^-'^' case. 

We can immediately use the MA^"'"^ result to solve Eq. 



(20 1, but only up to the n = 2 level: V2 = Q2V1. To solve 



the n — 1 equation exactly we will need to construct a 



set of recurrence relations involving Eqs. (40) and (41) 



using Dyson's identity. The resulting equations take a 
form similar to Eqs. (31 1-(36) if we define 



5fn,7n{qi) = /n,m('?l) " ^(1 - e"^')/n, 

+ 3(1 



-iqi 



n,m+l ; 



(43) 



for the 5fn m functions have precisely the same form 
as Eqs. (31l-(36l with n ^ n — 1 and cu ^ u! — il. 



As before, we write this set of equations in the form 
(5v„ = AJj(5v„_i +B;j(5v„+i, where 

6Vn = {Sfn,a, ■ ■ • , Sfn^n-l, 

Sfn.O^ ■ • ■ > ^fn,n-2^^fn.l' ■ ■ • : '^/rtn-l)^ > (47) 

and the solution can again be written in terms of an 
infinite continued fraction of matrices: 



5v„ = R„(5v„_i = 



1 



-A'JVn-l, (48) 



1 _ W A' 

" 1 _ . . . "+1 



where A[^{oj) = A„_i(tj — il) and Bjj(cj) = B„_i(a; — fi). 

To finish the calculation we will require explicit expres- 
sions for , fl^l , fiQ^^ , and fi}^^ . Using the definitions 



show that 



in Eqs. (27), (28), (40), and (41 1 it is straightforward to 



/M = ^EGo(fc-'Zi,--17)^-L_ 

91 



X [/2,o(gi) - /2.i('?i)] + eo - eie-^'=", (49) 



/S = -^EGo(fc~'Zi,--f^)^^^ 

91 



X [/2,o(gi) - /2,l('?l)] - eo + 



eie 



(50) 



^/n,o(<Zi) = /„,o('7i)-g(l-e''^)/„,o 
+ 5(1 - e-'«0/«,i - 5e-'('-^^^ni " e''')fn,n-i, (44) 

<5/n,n-l('?l) = .frun-l{qi) - .g(l - e'"^^)fn,n-l 

+ gil - e-"^0/n,n +5e'^'-^^'(l ~ e-^«0/+i' (45) 
for n > 1, and 

SfiA^i) = /i,o('Zi) -5(1 - e^'O/l'o^ +5(1 - e-^'^)/l\^ 

_g^~'^ik^,i)^l „ e»«0/ro'^ +5e'^'-«^)(l -e-^^O/i""!'^- 

(46) 

We have again added the superscript (1) to distin- 
guish /l,o\/l,\\/i:^'\/itl'^ from the MA(--°) expres- 
sions. With these definitions, the recurrence relations 



and 



91 

; [/2,o(gi) - /2,i('7i)] 



1 - e"?i 
el-e2e-^'=^ (51) 



91 

[/2,o(9i) - hAqi)] - ei 



1 



1 - e-^-Ji 



^626*'^°. (52) 



To combine all of these results and obtain a system of four 
equations in the four unknowns q' , /|\^ , q^-* , /^[^"^ , we 

need to rewrite f2,o{qi) — f2,i{qi) in terms of these four 
quantities. Using the definitions in Eqs. (44)-(46l " — 
Q2V1, and (5v2 = Ri(5vi, we find that 



V2 = 



10 



hM - hiiqi) = G^\k -q,uJ- n)Go{k - q, Q) {(i?oo - Rw) [Go(fc - qi,uJ - n)\g,, 

+5(1 - e^«OQ2Vi|o - 25 [1 - cos((?i)] Q2V1I1 + 5(1 - e-*'0Q2Vi|2 

+ge-*(fe-9i)(l _ e^9i)Q2Vi|4 + ge^'^-'^'^l - e-'«i)Q2Vi|5} , (53) 



where a; = w — 17 — g{RoQ — -Rio) ■ In this notation Rij are 
the matrix elements of Ri and Q2Vi|i is the i*^^ element 
of the product between the matrix Q2 and the vector vi, 

i.e., Q2Vi|i — Qio/i.x) + Qiifi]i + Qi2fi,'a^ + Qi'ifi,\^\ 
where Qij are the matrix elements of Q2. Inserting Eq. 
(53 1 into Eqs. (49 1- (52 1, and performing the required 



sums over momentum, we obtain the desired four equa- 
tions in four unknowns. The explicit expressions for the 
momentum averages that appear in these expressions are 
given in Appendix [X] Once this system has been solved 
we can easily compute the self-energy as 



{k,u:) 



9U^ 



(54) 



Again the MA*^"'^) expression is slightly more involved 
than the zeroth order approximation because we need to 
evaluate two continued fractions, but it is still numeri- 
cally trivial to evaluate. 



IV. RESULTS 

We present results for the ID breathing-mode Hamil- 
tonian using the variational MA approximation. Com- 
parisons are made to ED data^^ where possible. 



result in Fig. ^^c). Again, when plotted over the full pa- 
rameter range, the results are nearly indistinguishable. 
A look at the percent difference indicates that the rela- 
tive error is less than 2.5%. This is a truly remarkable 
fact considering that the qp weight contains information 
on the nature of the eigenstates, something that is rarely 
obtained accurately when using approximate methods. 

It is clear from Fig. |5]that the variational MA approx- 
imation cures all of the shortcomings of the simple MA 
generalization discussed in Sec. |IID[ A comparison of 
the variational MA and ED results shows that the MA 
energies are slightly higher than the exact numerical re- 
sults, as expected from a variational method, and that 
the asymptotic behavior predicted from both the pertur- 
bational theory and ED result is reproduced. These suc- 
cesses were expected because the variational MA approx- 
imation was designed precisely to remedy these problems, 
but it is still very encouraging to see that the physical 
picture described and used to motivate the variational 



MA approximation in Sec. Ill was indeed correct. 

In Fig. [6|^a) we plot the average number of phonons 



in the cloud, N„h, as calculated from the Hellmann- 
Feynman theorem.lIS] Since there is no numerical data 
available for this quantity we compare our findings to the 
standard Rayleigh-Schrodinger (RS) perturbation theory 



A. Ground State Properties 

In Fig. [5] we plot both the ground state energy and 
qp weight as a function of the electron-phonon coupling 
strength, using the variational MA approximation. The 
variational MA results show a clear improvement over 
the MA results shown in Fig. |2] and the agreement with 
the numerical data is excellent. The approximate and 
exact numerical data results are indistinguishable when 
plotted over the full parameter range of Fig. [IJa). To 
gain a closer look at the success of our approximation we 
plot the relative error between the MA results and the 
ED results in Fig. |5jb). Indeed the approximation is 
giving extremely close agreement with the numerical re- 
sults, with less than 0.3% relative error for both MA^^'^^ 
and MA^"'!). The largest errors occur at intermediate 
couplings, as expected because the MA approximation 
is exact in both the zero coupling and zero bandwidth 
limits. We also show a comparison of the quasiparti- 
cle weight calculated using variational MA with the ED 
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A.=g7(nt) 



0.5 1 1.5 

A.=g7(nt) 



FIG. 5: (Color online) (a) Ground state energy, (b) Percent 
difference from ED ground state energy results, (c) Ground 
state qp weight, and (d) Percent difference from ED qp 
weight results, as a function of the effective coupling A, for 
t = 1, Q./t = 0.5. The ED results are from Ref. [T71 
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again very good, coming well within 0.5% relative error 
of the numerical result, showing that the variational MA 
approximation is accurate for all momenta k. The rea- 
son for the non-monotonic polaron dispersion has been 
discussed at length in Ref. XL ^^^d is due to a larger 
effective 2nd nearest neighbor hopping than the effective 
nearest neighbor hopping of the polaron - this is a direct 
consequence of the structure of the polaronic cloud. In 
Figs. [tJ^c) and |7|^d) we plot the qp weight and its rela- 
tive error as a function of momenta. We again find good 
agreement. As explaned, the agreement is expected to 
improve for both smaller and larger A values, where our 
approximations become asymptotically exact. 



FIG. 6: (Color online) (a) Average number of phonons and 
(b) effective mass, (c) Percent difference from ED effective 
mass results, as a function of the effective coupling A, for 
t = 1, n/t = 0.5. The ED results are from Ref. HIl 



at small couplings, and to the strong-coupling pertur- 
bation theory for larger couplings. The approximation 
reproduces both asymptotic limits to a high degree of 
accuracy, as should be expected based on the success of 
the ground state energies and qp weights shown in Fig. 
[5] Since we are unaware of any explicit expression for 
the RS energy for the breathing-mode Hamiltonian in 
the literature, we state the result here. By evaluating 
the standard expression for the first order correction 



E, 



(1) 



(55) 



one can easily show that 



(1) 



2.9 5o(f^-efc)- Y cosfc [1 -{n- ek)goi^ - ^k)] 

(56) 

for the ID breathing-mode Hamiltonian. 

In Fig. |6]we also plot the effective mass as a function 
of the coupling strength. Again the agreement with the 
numerical data is excellent, as confirmed by the relative 
errors shown in Fig. [6jc). 

B. The Polaron Band 

With our analytical expression for the self-energy we 
can also calculate momentum-dependent results. In Fig. 
|7|^a) we plot the lowest energy state for momenta < 
k < n and compare our results to the available numerical 
data. We again find that the variational MA approxima- 
tion is highly accurate. Because the polaron dispersion 
is relatively flat for the intermediate coupling strength of 
A = 1.071502 shown here, the energy range of the plot is 
quite narrow and we can clearly discern the difference be- 
tween the approximate and numerical results. However, 
by looking at the relative error in Fig. [7|b), we demon- 
strate that the accuracy of the variational MA result is 



C. Higher Energy Properties 

Lastly, we consider the high energy properties of 
the ID breathing-mode Hamiltonian using and 
MA(^'i). In Fig. H we compare our predicted spectral 
weights to available numerical data. When plotted on a 
linear axis, the results are essentially indistinguishable, 
especially for MA^^'^). To gain a better view we dis- 
play the same plots on a logarithmic scale in Fig. [sf^b). 
There are a few interesting features to note. At the 
MA^^^ level (red dashed line) the majority of the spec- 
tral weight is found to be in the correct location, however, 
the electron-l-phonon continuum is completely absent for 
the parameters shown, as expected. At the MA^"'^^ level 
the continuum is reproduced in the expected location, 
in very good agreement with the ED prediction. Fur- 
thermore, the finite size effects responsible for the sharp 
peaks in the continuum of the ED result are absent from 
the MA^"'^) data. As a guide illustrating the correct lo- 
cation for the polaron+phonon continuum we have added 
the vertical blue dashed lines to denote Egs and Egg + U 
for the MA^'"'-'^) ground state energy. The lower edge is 
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FIG. 7: (Color online) (a) Polaron dispersion Ek, (b) Percent 
difference from ED ground state energy results, (c) quasi- 
particle weight Zk, and (d) Percent difference from ED qp 
weight results. Results are shown for t = l,Q/t — 0.5, and 
A = 1.071502. The ED results are from Ref. El 
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FIG. 8: (Color online) (a) A{k = 0,uj) vs. lo and (b) lnyl(fe = 
0, io) iort = l,n = 0.5, T] = 0.004, and A = 0.75, 1.0, 1.25, 1.5; 
The vertical blue dashed lines denote Egs and Egs + using 
MA'"'^'. The ED results are from Ref. El 
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FIG. 9: (Color online) (a) A{k,uj) and (b) lnA(k,uj) vs. ui for 
fc/TT = 0, 0.25, 0.5, 0.75, 1, t = 1, = 0.5, A = 1.25, ri = 0.004; 
The ED results are from Ref. [13 



located correctly, however especially for smaller A, MA 
predicts a somewhat wider continuum than ED. Given 
the very limited availability of numerical results of this 
type, we do not know if this discrepancy is due to trun- 
cation approximations in the ED solution, or is due to 
inaccuracies of the MA approximations. 

In Fig. [9] we plot the spectral weight for a fixed cou- 
pling strength and vary the momentum k. Again on the 
linear scale the results are very encouraging. The MA 
result, particularly MA'"'^) predicts the correct location 
for the spectral weight over the entire energy range. How- 
ever, a logarithmic scale plot does reveal a notable short- 
coming of the approximation. For fc = the agreement 
with the numerical data is excellent and the continuum 
is located at the expected Egs + fi. For k — 0.257r and 
k = O.Stt the agreement is still good, but as we begin 
to approach the band edge we see a significant deviation 
of the MA result from the ED result, with the MA con- 
tinuum coming to much too low energies. This feature 



has extremely little spectral weight (see Fig. [9|a)) and is 
extremely unlikely to be detectable in any experimental 
realization, however it is in stark disagreement with the 
expected result confirmed by the ED data. We do not 
currently understand this failure at higher k values. 

Finally, in Fig. [TO] we show more detailed plots of 
the evolution of the spectral weight for fc = from zero 
coupling to large coupling, for a full range of energies. 
The two panels compare the MA^"-") (a) and MA^''-^) 
(b) solutions. The one obvious difference is the location 
of the polaron-l-phonon continuum, which is incorectly 
located around —2t for MA^'"'^\ whereas it always starts 
at fl above Egg for MA''"^^^ Aside from this feature 
with rather little weight, all the features which contribute 
significantly to the spectral weight are in good agreement. 

Plots of this nature would require extensive numerical 
computational time if one tried to generate them using 
"exact" computational methods. For example, achieving 
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FIG. 10: (Color online) Spectral weight A{0, uj) vs. uj in ID 
using (a) MA*'''^' and (b) MA*"'^'. The resuhs are shown for 
t = l,n — 0.5, ri = 0.01, and A varying from to 2. Curves 
corresponding to A — 0.3,0.6,1.2 and 1.8 are highlighted in 
red. 



convergence for the ED results at larger A required inclu- 
sion of billions of states, in some cases, in the truncated 
Hilbert space. Diagonalizing such large (even though 
sparse) matrices is not a trivial task. By comparison, be- 
cause we have an analytical expression for the self-energy, 
these detailed MA plots covering many energies and cou- 
plings take just seconds, at most minutes, to generate. In 
fact, MA is much faster even than SCBA, which for this 
model requires many numerical integrals over the Bril- 
louin zone (and is of questionable use for medium and 
large couplings). 



V. SUMMARY AND CONCLUSIONS 

In summary, we have presented a way to general- 
ize the momentum average approximation to more gen- 
eral electron-phonon coupling models with momentum- 
dependent couplings. The approximation still sums all 
of the self-energy diagrams, albeit with each diagram ap- 
proximated in such a way that the full sum can be per- 
formed, and it is exact in both the zero coupling and 
zero bandwidth limits. As in the application of MA to 
the Holstein model, the approximation is analytical and 
easy-to-use, and gives highly accurate results over the 
entire parameter space. In this paper we have actually 
presented two different generalizations of MA, a straight- 
forward generalization that can easily be applied to any 
model with a momentum-dependent coupling with mini- 
mal effort, and a more specific (and extremely accurate) 
generalization that takes the details of a given model into 
account and enlarges the variational subspace accord- 
ingly. We have used the ID breathing-mode Hamiltonian 
as an example to gauge the accuracy of both of these MA 
generalizations. While the straightforward extension of 
MA gave fairly accurate results, it was non-variational 
and it approached the zero-bandwidth asymptotic limit 
very slowly. The variational MA approximation remedied 
both of these problems and produced extremely accu- 
rate results, coming well within 0.3% error of the avail- 
able numerical results for the ground state energies of 
the breathing-mode polaron. We also showed that MA 
can be systematically improved in both cases, leading 
to higher accuracy and the correct location for the elec- 
tron+phonon continuum. 

The successful generalization of MA to this much 
broader class of models in very encouraging. Numer- 
ical studies of models as complicated as even the ID 
breathing-mode Hamiltonian are very intensive, and the 
MA approximation provides a quick and easy-to-use way 
to gain an understanding of these more realistic models 
without having to do a detailed numerical analysis from 
the start. We hope that this tool will be extremely useful 
for probing more realistic and experimental realizations 
of interesting physical systems. The range of applicabil- 
ity of the MA approximation has been growing steadily. 
As mentioned previously, it has been successfully applied 
to systems with multiple phonon modes and multiple 
free-electron bands. Other obvious and very useful gen- 
eralizations of this approximation would be to apply it to 
finite particle densities and/or finite temperatures, and 
possibly to an even broader class of Hamiltonians involv- 
ing electron-electron interactions. 
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breathing-mode model one can show that 



APPENDIX A: MOMENTUM AVERAGES 



go{k,uj) = 2g goM 
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In this section we derive the various momentum av- 
erages required to evaluate the self-energy expressions 
derived in this work. We require the following weighted 
momentum averages of the free electron propagator for a 
ID tight-binding dispersion: 
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It is straightforward to show that 
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In the MA^""'!) calculation we require the momentum 
average of Eq. ( 53 1 . For the ID breathing- mode Hamil- 



tonian (i.e. gq — — 2i(7sin(g/2) and Sk = — 2tcos(fc)) we 
require the following momentum averages: 
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ff2(w) = -goi^j ^ — 9i{^j, (A4j and 
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For the MA^^^ and MA'^^ self-energy expressions we 
also require 5o(k, w), as defined in Eq. (15 1. For the ID 
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